NEW, HIGHLY ACCURATE PROPAGATOR FOR THE 
LINEAR AND NONLINEAR SCHRODINGER EQUATION. 



HILLEL TAL-EZER*, RONNIE KOSLOFFt, AND IDO SCHAEFER * 

Abstract. A propagation method for the time dependent Schrodinger equation 
was studied leading to a general scheme of solving ode type equations. Standard space 
discretization of time-dependent pde's usually results in system of ode's of the form 



where G is a operator ( matrix ) and m is a time-dependent solution vector. Highly accu- 
rate methods, based on polynomial approximation of a modified exponential evolution 
operator, had been developed already for this type of problems where G is a linear, time 
independent matrix and s is a constant vector. In this paper we will describe a new al- 
gorithm for the more general case where s is a time-dependent r.h.s vector. An iterative 
version of the new algorithm can be applied to the general case where G depends on t 
or u. Numerical results for Schrodinger equation with time-dependent potential and to 
non-linear Schrodinger equation will be presented. 
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1. Introduction. The time dependent Schrodinger equation is of fun- 
damental importance, it governs quantum dynamics. As a result any sim- 
ulation of quantum phenomena requires an effective scheme to represent 
and solve this equation: 



where V' is a vector representing the wave function and H the Hamiltonian 
operator [20] • Applications differ considerably. The dimension of Hilbert 
space required to represent the wave function ^ can vary from 2, for a two- 
level-system, to ~ 2'^" in practical applications. When the size of Hilbert 
space becomes too large to be represented directly, approximate methods 
are employed which lead to a nonlinear version of the Schrodinger equation 



The central role of Eq. (jl.ip in quantum dynamical simulations has 
generated a wealth of numerical methods to solve the equation. For low 
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dimensions, the common approach is based on diagonahzing the Hamil- 
tonian operator H . For higher dimensions, this becomes impractical and 
one has to resort to matrix-free methods which require only the evaluation 
of the operation of the Hamiltonian on a vector. As a result, an implicit 
knowledge of the Hamiltonian is sufficient. 

Many methods have been developed and implemented to propagate 
the equation in time. Typically the propagation period is divided into time 
steps. As a result the error in each time step will accumulate. This means 
that effective methods should have as large as possible time step and have 
a high accuracy within a time step. For time independent Hamiltonian 
operators a global polynomial expansion of the propagator is the method 
of choice [3]: 

(1.2) i^it) = e-'^V(O) - J2 ^n{t)Pn{H)^m 

n 

where P„ (x) is a polynomial of order n which is evaluated recursively. The 
most popular choice has been the Chebychev polynomial [T^ due to its 
exponential rate of convergence. Other polynomials have been tried with 
similar or inferior results. 

In many applications, the Hamiltonian is explicitly time dependent. 
These include systems subject to a time dependent electromagnetic field 
(spectroscopy), quantum control which requires to infer the time dependent 
field that leads to a desired outcome such as a quantum gate [7]. In these 
problems the remedy to overcome the explicit time dependence was to 
employ a short time step in which the field is approximated as piecewise 
constant. This solution immediately degrades the accuracy to first order 
in the time step. Four general approached have been explored to overcome 
this difficulty. 

1. Solving the equations using general Taylor based solvers such as 
Runge Kutta or second order differencing. These methods have 
slow convergence properties [8]. 

2. Employing the {t,t') method which eliminates the explicit time 
dependence by embedding the problem in a larger Hilbert space 
adding time translation to the Hamiltonian H{t) — > H{t') + i-^- 
The method restores the accuracy of the high order polynomial 
expansion but has been found to be expensive in use[T^. 

3. Another class of approaches rely on the Magnus expansion to over- 
come the problem of time ordering [T^. The solution is cast into the 
form: tp{t) = e^ipi^) and approximated as « e"^^ e^^ e^^ . 
This type of solution includes the split operator method[T5] as well 
as polynomial approximations of the exponent [1]. 

4. When the Hamiltonian can be split as: H = Hq + V{t) then 
iP{t) = e-*^°*V'(0) - i e-~'"'>^'-*'^Vit')iljit')dt' . This formal so- 
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lution establishes the base for a polynomial approximation of the 
result [in], SI]- 

When considering nonlinear version of the Schrodinger equation, such 
as the GrossPitaevskii equation or time dependent density functional equa- 
tions, methods 2 and 3 are not applicable and we are left with options 
based on 1 and 4. The new algorithm presented in this paper belongs 
to the fourth approach. We will demonstrate that the new algorithm is 
highly efficient with respect to accuracy versus numerical effort, both for 
linear time dependent problems as well as for non linear versions of the 
Schrodinger equation. 

2. The new algorithm (linear case). Let us consider a general 
system of ode's of the form 



where G is a constant, N x N matrix. If s is constant then, by Duhamel 
principle, the solution is 



(2.1) 



ut = Gu + s. 



(2.2) 



w(0) = Wo, 



(2.3) 




Formal integration results in 



(2.4) 



M(t) = e^*«o + /i(G,t)s, 



where 



(2.5) 




Since 



(2.6) 



e 



,zt 



zh{z,t) + l, 



then 



(2.7) 



G/i(G,t) + / 



and therefore 



(2.8) u{t)^vo + h{G,t)v^ 
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where vi = Gvq + s. 

Going one step further, let us consider the system 

(2.9) ut = Gu + So + tsi. 

Using similar steps as above, wc get the formal solution 

(2.10) u{t) =vo+ tvi + /2(G, t)v2, 
where 



(2.11) f2{z,t) = 



1 l„zt 



(e^* -l-zt) ^ ^ 
^ . = 



and V2 = Gvi + si. The following Lemma applies to the general case. 
Lemma: 

The formal solution of the set of ode's 

m-l • 

(2.12) „^ = Gn+^-s,- 
is 

m-l 

(2.13) U=J2 + fm {G, t) Vm, 

where vj satisfy the recurrence relation 

(2.14) vo = uo 

(2.15) Vj = Gvj-i + Sj-i ^ < j < m 
and 

(2.16) fm{z,t) = \ ^-l' ^^=0 i^- ) 
Proof: 



It is easily verified that 
(2.17) % = ^/" ' 



dt •"" (w-l)! 
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Hence 



(2.18) = 51 ~T^J + l + Gfr,^V„, + 



V. (m — 1)! 

or 

m-l 

(2.19) = ^ — + 1 + GfrnVm ■ 

Using (|2.15p we get 

m — 1 m — 1 

(2.20) ut = + ^ -s, + G/,„w„ 

3=0 ^' j=0 ^' 

Hence 

^m-l \ m-l 



(2.21) Ut = G ^ -fj + fmVrn + J] " 



or 



J- ]■ 



?7l— 1 



(2.22) = Gm + ^ — 

and the proof is concluded. 

remark: When z is very small, computing (z,t) as defined in ()2.16p 
can be unstable due to roundoff errors. Possible remedy is to use instead 
an approximation based on Taylor expansion 



(2.23) Uiz,t) = t"^J2 



[zty 



The solution vector u can be approximated with high accuracy as 

m-l 

(2.24) w » ~! + ^"^ 

where pk{z,t) is 'optimal' polynomial which approximates fjn{z,t) where 
z G D and D is a, domain in the complex plane which includes all the 
eigenvalues of G. The pk polynomial can be based on Chebyshev expansion 
[TBI , Arnoldi approach [4l, [T7l or Newton interpolation approach ^ISj . 
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In the more general case where s is any function of t we do first Cheby- 
shev approximation of s 

m — 1 

(2.25) s(t)« ^S,T,(t) 

and then transform the expansion to the Taylor-like representation as in 

3. Time-Dependant G. Let us consider now the case where the ma- 
trix G depends on t 



(3.1) Ut^G{t)u + s{t), u° = u{0), 0<t<T. 

(The time dependent Schrodinger equation where the potential depends on 
t is an example of such an equation). 

In order to apply the new algorithm in this case, one has to resort to a 
time-steps algorithm. Consider that we have marched already to time level 
tn and we want to compute the solution at time level tn+i- p.ip can be 
written as 

(3.2) Ut = GnU + Snit), 0<t<T 

where 

(3.3) G„ = G f i„ + — j , s„ (0 = s{t) + {G{t) -Gn)u, 
and 

(3.4) At = i„+i - i„. 

Observe that s„ (t) depends on u which is unknown yet at the time interval 
[tn,tn + At] but, as described in Main Algorithm below, a set of approx- 
imated vectors w:^ , j — 1 , . . . , m which approximate the solution at the 
Chebyshev time points 

(3.5) tj =tn + ^{l- Vj) , Vj = cos (^^^^^^3Y") ' 1 < « < "I, 

can be computed in the previous time step and is used to compute the 
Sj vectors as defined in (|2.12p . Only in the first step one has to use an 
iterative algorithm in order to compute the set of approximated solution 
vectors at the time points 

(3.6) t, = —{l-y^), y, = cos ' , l<J<m 
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where the first guess is 

(3.7) u] = 7/°, 1 < j < m. 

The iterative algorithm is stopped when — w^J| satisfies the desired 

accuracy. 

First Step Algorithm 

Given: u°, e, m and let tj ^ ^ (l - cos , l<j<m 

1. ) Uj =u", j = l,...,m 

2. ) Compute Sj = so{tj) (defined in (13.31) ) 

3. ) Compute Sj (defined in (I2.12p ) by using cosine transform of Sj and 

then Taylor-like transform 

4. ) Use (pHl) to compute uf^"^, l<j<m 

5. ) if ~ ""mil — ^ then stop 

6. ) Uj = uf^"^, l<j<m 

7. ) go to 2 



After computing the initial solution vectors at the time points tj = 
^1 — cos ^-^"^^^^^ ^ , 1 < J < we are ready to continue with the main 
algorithm which computes the solution at the time interval [0,r]. 



At 
2 



Main Algorithm 

Given: vq, m, {-Ujjjii, T 
t = 0,n = 

1.) Let t] ^ t+M (i _ cos ((i^)) , = t+^t+M (i _ cos (%^)) 

1 < J < TO 

1. ) Compute Sj — s„(tj), (defined in p.Sp 'l 

2. ) Compute Sj (defined in (I2.12p ) by using cosine transform of Sj and 

then Taylor-like transform 

3. ) Use (i2T4| to compute at {t]}]Li 

4. ) \it = T then stop 

5. ) t = t + At, + 1 
7.) go to 1 

Observe that t}-^^ = tl = t + At hence, at each step, the solution vector 
at this point is computed twice. The first one is the predictor and the 
second one is the corrector. 

4. Nonlinearity. Let us consider now the nonlinear case. 
(4.1) ut = G{u)u + s{u), 0<t<T. 
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Implementation of the new algorithm in this case is almost the same as it 
is done in the case described in the previous section. (j4.ip can be written 
as 

(4.2) Mt = G„u + s„, 0<t<T 
where 

(4.3) G„ ^g(u (^t„ + y s„ = s{u) + {G{u) - G„) u. 

The rest of the description of the algorithm is exactly the same as in the 
previous section. 

5. Numerical Examples. The numerical examples presented in this 
section address the case where the eigenvalues of the spatial matrix G are 
on the imaginery axis. In this case one can use the Chebyshev approach 
[16] . In a future paper we will treat the more general case ( e.g. boundary 
value problems, advection diffusion ) where the domain of eigenvalues is on 
the left side of the complex plane. 

Example 1: Time-dependent r.h.s 

Let us consider the differential equation 

(5.1) ut = Ux + s{x, t) < a; < 27r 

where 



s{x, t) ~ sin(6a;) cos(t)— 2 cos(lOx) cos(2t)— 6 cos(6a;) sin(t)+10 sin(10a;) sin(2t). 
(5.2) 

The exact solution is 

(5.3) u{x, t) = am{t) sin(6a;) + sin(2t) cos(lOx). 

Since we have periodicity in space we can use spectral Fourier for space 
approximation. It results in a set of ode's 

(5.4) ut^Gu + s 

where u is a vector of length n ( number of grid points ) , G is an n x n matrix 
which carries out the Fourier spectral differentiation and s is a vector of 
length n which is time-dependent. 

We have solved this problem in the time interval [0 5]. Since the 
solution is periodic with highest mode equal to 10, using n = 32 is suffice 
to compute exactly the spatial derivative. Hence, the error comes solely 
from time approximation. 
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In order to compute solution in this time interval which satisfies 

(5.5) ||u - UexactW < 10"^ 

we had to use m = fc = 14 ( these parameters are defined in (I2.12[) 
and (|2.14l) ). It means that all together we had to do 28 matrix- vector 
multiplications. 

Applying standard ODE45 for this problem, we had to do 860 matrix- 
vector multiplications in order to compute the solution to the desired ac- 
curacy. 

Example 2: Schrodinger equation with time-dependent potential 

As a second example, we consider harmonic oscillator of mass m ~ I 
and frequency lu — 1 driven by a linearly polarized electromagnetic field 
with frequency v — 1. We have to solve 

(5.6) ipt^-iH{r,t)^P 
where the time-dependent Hamiltonian is given by 

(5.7) Hir,t) = -i^ + ir^ -f rsin^ (^^) cos(t) . 

The final time is set to T = 15 . The Hamiltonian is represented on a 
Fourier grid with n = 128 grid points, and Tmax = 10 = —rmin- We have 
used the spectral Fourier method to approximate the spatial derivatives. 
Taking the initial wave function to be 

(5.8) V(»',0) = e-"^' 

we computed the numerical solution by two methods : 

1. ) RK4 (Runge-Kutta of order 4 ) 

2. ) the new algorithm. 

In all the tables below, matvecs represents the number of applications 
of the Hamiltonian. 

The first table presents the RK4 results. For stability , the time step 
should be At = 0.01, hence the minimal number of time steps needed to 
march to T = 15 is 1500. 
Table 1-RK4 



time-steps 


matvecs 


Relative L2 Error 


1500 


6000 


5.6e-04 


3000 


12000 


3.5e-05 


6000 


24000 


2.2e-06 
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Observe that dividing the time step by 2, the error is reduced by a 
factor of almost 16 as it should be since RK4 is a scheme of order 4. Hence, 
in order to get high accuracy, e.g. of order 10"^" , one should do around 
192000 matrix-vector multiplications. 

In the next few tables we present the results for the new algorithm. 
The tables differ by the m and k parameters where m is the number of 
Chebyshev points in the interval [t,t + At] and k is the degree of the 
polynomials used to approximate the function /„. 
Table 2- new algorithm, m=k=7 



time-steps 


matvecs 


Relative L2 Error 


350 


4563 


3.7e-02 


400 


5213 


3.9e-08 


600 


7813 


3.9e-10 



Table 3- new algorithm, m=k=8 



time-steps 


matvecs 


Relative L2 Error 


300 


4515 


2.3e-08 


400 


6015 


1.3C-09 


450 


6765 


4.8e-10 



Table 4- new algorithm, m=k=9 



time-steps 


matvecs 


Relative L2 Error 


280 


4777 


6.0e-08 


350 


5967 


1.4C-09 


400 


6817 


3.2e-10 



Remark: The minimal number of time-steps presented in the last 3 
tables were such that taking smaller number will result in instability. 

Observe that the new algorithm is significantly more efficient then RK4, 
especially when one is interested in high accuracy. In this case, the new 
algorithm is almost 30 times more efficient then RK4. 

Example 3: Nonlinear Schrodinger equation 

For a nonlinear example we choose the GrossPitaevskii equation de- 
scribing the dynamics of a Bose-Einstein-Condensate (BEG) in a harmonic 
trap: 



(5.9) 



—iH (r, tp) ip 
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where the Hamiltonian is given by 



(5.10) H{r,^) = -l^ + \r' + \^^\\ 

The final time is set to T = 10. The Hamiltonian is represented on 
a Fourier grid with n = 128 grid points, and fmax = ^y/n = — rmin- The 
spectral Fourier method is used to approximate the spatial derivatives. 

The initial state is 

(5.11) iJo = e'^'^vo 

where vq is the eigenvector related to the smallest eigenvalue of the non- 
linear Hamiltonian. 

As in the previous example, we computed the numerical solution by 
RK4 and by the new algorithm. 

The next table presents the RK4 results. For stability , the time step 
should be At = 0.01515, hence the minimal time steps needed to march to 
T = 10 is 660. 

Table 5-RK4 



time-steps 


matvecs 


Relative L2 Error 


660 


2640 


4.96e-01 


1320 


5280 


2.00C-02 


2640 


10560 


1.20e-03 


5280 


21120 


7.29e-05 



Taking into account that RK4 is a scheme of order 4 we can conclude 
that in order to get high accuracy, e.g. of order 10~^° , one should do 
around 382000 matrix-vector multiplications. 

In the next few tables we present the results for the new algorithm. As 

in the previous example, the tables differ by the m and k parameters. 

Table 6 - new algorithm, m=k=7 



time-steps 


matvecs 


Relative L2 Error 


300 


4043 


3.3e-05 


500 


6643 


9.5e-08 


700 


9243 


1.7e-09 
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Table 7- new algorithm, m=k=9 



time-steps 


matvecs 


Relative L2 Error 


200 


3587 


5.67e-05 


300 


5287 


1.14e-07 


400 


6987 


3.00e-09 


500 


8687 


6.43e-10 



Observe that for moderate accuracy of order 10^^, 21120 matvecs were 
needed in the RK4 case while using the new algorithm with m = k = 9, only 
3587 matvecs were needed. The increase in efficiency is more pronounced 
when high accuracy is needed. For order of 10^^" accuracy, 382000 matvecs 
are needed in the RK4 case compared to 8687 matvecs for the new algo- 
rithm. 

6. Conclusions. In this paper we have presented a new algorithm for 
solving a class of linear and nonlinear Schrodinger equations which can be 
applied to general system of ode's. In the stationary linear case it is possible 
to reach the upper time level in one step with very high accuracy. Due to 
the fact that there is only one step, the accuracy is not deteriorating since 
there is no accumulation of errors. In the case where the matrix involved 
depends on time or in the case of nonlinearity, the time interval should be 
divided to time steps but the size of the time step is significantly larger 
than what is needed in standard explicit algorithms like Runge-Kutta. 

The high accuracy (spectral) of the algorithm can be traced to the fact 
that the algorithm does not use any Taylor considerations. Taylor theorem 
is an extremely important tool in analysis but due to its locality it can 
lead to inferior numerical approximation. We believe that whenever it is 
possible to develop an algorithm which is Taylor free, one should explore 
this possibility. 

7. Remarks. During the refereeing process we came to know of meth- 
ods known as Exponential Integrators (e.g. [5], [5], [3] ) which also make 
use of the functions defined in 12.161 The algorithms described in those pa- 
per are, like Runge-Kutta approach, based on Taylor considerations while 
the algorithm described here is Taylor-free. This is the main difference 
between the two approaches. Since the present algorithm is Taylor-free, 
there is no meaning to the term - "order of the method" which we have in 
the Exponential Integrator methods. 
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